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ABSTRACT 

We calculate the flux of internal gravity waves (IGWs) generated by turbulent convection 
in stars. We solve for the IGW eigenfunctions analytically near the radiative-convective in- 
terface in a local, Boussinesq, and cartesian domain. We consider both discontinuous and 
smooth transitions between the radiative and convective regions and derive Green's functions 
to solve for the IGWs in the radiative region. We find that if the radiative-convective transi- 
tion is smooth, the IGW flux ~ F^oaAdl H), where Fconv is the flux carried by the convective 
motions, d is the width of the transition region, and H is the pressure scale height. This can be 
much larger than the standard result in the literature for a discontinuous radiative-convective 
transition, which gives a wave flux ~ /^convAl, where M is the convective Mach number. 
However, in the smooth transition case, the most efficiently excited perturbations will break 
immediately when they enter the radiative region. The flux of IGWs which do not break and 
are able to propagate in the radiative region is ~ F^om-Mr'^^idl H)^^^ , larger than the discon- 
tinuous transition result by {MHIdY^^^. The transition region in the Sun is smooth for the 
energy-bearing waves; as a result, we predict that the IGW flux is about five times larger than 
previous estimates. We discuss the implications of our results for several astrophysical appli- 
cations, including IGW driven mass loss and the detectability of convectively excited IGWs 
in main sequence stars. 
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1 INTRODUCTION 

Internal gravity waves (hereafter, IGWs) are a class of waves 
in a stably stratified background in which buoyancy serves as a 
restoring force. IGWs propagate in radiative zones in stars and 
can influence composition, angular momentum, and energy trans- 
port within stars. IGWs could also be important diagnostics of 
stellar structure — the detection of standing I GWs fe-modes) has 
been a long-standing goal of helioseismology ( S evernvi et ahlig?^ : 
ISrookes et alj|l97q) . as g-modes provide better information about 
the core of the Sun th an the more easily obser ved global sound 
waves (p-modes) (e.g.. lTurck-Chieze et alj|200lh . However, IGWs 
are evanescent in the convection zone, so their surface manifesta- 
tion is expected to be small. 

IGWs have been invoked to explain the observation that 
F- stars have a smaller t han expected Li abundance (e.g., 
IXalon & Charbonnell 1 19981) . iGarcia Lopez & Spruitl ( Il99lh . here- 
after GLS91, first suggested that mixing from IGWs could 
enhance diffusion of Li, leading to lower Li abundances. 
ICharbonnel & Talorj ( |2005|) invoke IGWs to explain both the Li 
abundances of solar-type stars and the rotation of the solar interior. 
When propagating through a differentially rotating star, selective 
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damping of modes can deposit the wave 's angular momentum and 
modify the star's rotation profile (e.g.. I Kumar & Ouataerj [79971 : 
IZahn et alJl997l : lTalon et al.l2002l) . Note, however, that IGWs gen- 
erally have an anti-diffusive effect, accentuating angular velocity 
gradients. This anti-diffusive behavior leads to the quasi-biennial 
oscillation (QBO) in the Earth's atmosphere, and ha s been studied 
extensively by the atmospheric science community ( iBaldwin et alj 
l200ll : lFritts & Alexandeij|2003h . 

Massive stars have convective cor es surrounded by a radia- 
tive envelope. lOuataert & Shiodd j2012l) suggested that extremely 
vigorous convection within the last ~ year of a massive star's life 
could generate a super-Eddington IGW flux and drive significant 
mass loss. Earlier in a massive star's life, the angular momentum 
carried by IGWs may generate substantial differential rotation, per- 
haps mirroring the QBO in the Earth's atmosphere dRogers et alj 
I201I) . 

In some stars , IGWs are linearl y unstable, driven by, e.g., the e 
or K mechanisms jUnnoetalJI 19891) . Even absent such linear driv- 
ing, however, IGWs are thought to be generated by turbulent con- 
vection. Although IGWs are evanescent in a convective region, they 
can be excited by Reynolds stresses or entropy fluctuations associ- 
ated with the convection. A related excitation mechanism is gener- 
ation by overshooting convective plumes which penetrate into the 
radiative region. Numerical simulations of a radiative zone adja- 
cent to a convection zone find efficient generation of IGWs (e.g.. 
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Rogers & Glatzmaieil l2005h iMeakin & Amed 120071: iBrun et al 
201 ih . Al though simulatio n s rep orted in [Rogers & Glatzmaiei 
1 2005h & iMeakin & Arnetll ( l2007h show power distributed over 
a wide range of freq uencies and wavelengths, the power spectra 
m iBrun et alj |20l 1) exhibit ridges corresponding to discrete g- 
modesQ Simulations require artificially high diffusivities in the ra- 
diative zone to maintain a strong convective flux, and thus IGWs 
are artificially strongly damped in the radiative zone. This makes it 
difficult to estimate the IGW flux or quantitatively study the effects 
of IGWs on the stellar structure. 

There have been several efforts to analytically estimate the 
flux of IGWs stochastically excited by turbulent convection. These 
models are essential for determining the resulting efficiency of the 
mixin g, angular m omentum transport, or mass-loss produced by 
IGWs. lPressI jl98ll. hereafter P81) and GLS91 match pressure per- 
turbations in the convective region to pressu re perturbations in the 
waves, whereas iGoldreich & Kumaj 1199(1 hereafter GK90) and 
B09 calculate eigenmodes and derive how their amplitudes change 
using an inhomogeneous wave equation. P81, GLS91, and GK90 
all model the convective region using mixing length theory, assum- 
ing a Kolmogorov turbulence spectrum. B09 uses an energy spec- 
trum calculated from a direct numerical simulation of the solar con- 
vection zone. Each of these papers predicts a different IGW power 
spectrum. 

In this paper, we calculate the IGW flux generated by turbu- 
lent convection and clarify the relationship between different pre- 
dictions in the literature. In Section [2] we state our assumptions 
regarding the background state, and describe some properties of 
IGWs. Our main calculation is in Section [3] where we introduce 
our formalism for calculating the IGW flux. Our formalism relies 
on calculating a Green's function using the eigenmodes of the sys- 
tem (also discussed in FBI). We relate our method to GK90's in 
Appendix [C] In Section 13.51 we calculate the IGW flux and rms 
wave displacements for both smooth and discontinuous radiative- 
convective transitions. Next, we show that our results for a discon- 
tinuous transition can be derived more heuristically using pressure 
balance arguments (Section^; we also make detailed comparisons 
to previous results (Section|5}. Finally, in Section|6]we conclude, 
show how our results increase the predicted IGW flux in stars, and 
discuss some implications of this increased wave flux. 



2 BACKGROUND STATE AND PERTURBATION 
EQUATIONS 

In this paper we consider a simple model of a radiative zone and 
adjacent to a convection zone. We assume that the length scales of 
interest are small in comparison to the stellar radius, i.e., we are 
in the local limit, so we use cartesian geometry, where e, is the 
direction of gravity. In our model, the radiative zone is the region 
—L < z < Zi, and the convection zone is the region Zi < z < L, 
where z, is the location of the radiative-convective interface, and 
both regions have a horizontal area y{. We take L and to be 
much larger than any other length scale in the problem, and will 
assume z, is close to zero. In Figure [T] we sketch a schematic of 
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solve the anelastic equations, which do not conserve energy jBrownetal] 
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Figure 1. A schematic of our problem setup. The radiative-convective in- 
terface is at z = Zi, where Zi is close to zero, and has width d. Gravity points 
downward in the z direction. The convection zone is the region z > Zi and 
the radiative zone is the region z < Zi- We will use f^jad to denote the part of 
the vertical displacement within the radiative zone. If d is small, the waves 
see the radiative-convective transition as discontinuous; we will use super- 
script d to denote results for a discontinuous transition. If d is large, the 
waves see the radiative-convective transition as smooth; we will use super- 
script .s to denote results for a smooth transition. Eans.l4n& 1491 give our 
IGW flux estimates in the discontinuous and smooth limits, respectively. 



our model. Using a domain with finite vertical extent provides sim- 
pler boundary conditions, but yields the same results as an infinite 
domain. 

Furthermore, we employ the Boussinesq approximation. This 
is appropriate if the wave generation occurs close to the radiative- 
convective boundary, and if we are only concerned with IGWs near 
this boundary. We will see that the wave generation primarily oc- 
curs in a region with height approximately equal to the size of 
the energy bearing convective motions, which we assume is ~ H 
the pressure scale height. Although the Boussinesq approximation 
is only rigorously valid on length scales smaller than H, we re- 
cover results similar to those presented in GK90 who used the fully 
compressible equations. We thus believe that our results would not 
change qualitatively if we used the fully compressible equations. 

We model the radiative region as a stably stratified atmosphere 
with a squared buoyancy frequency N^. The convective region is 
much more complicated due to turbulent motions. We decompose 
the fluid properties in the convection zone into time averaged and 
fluctuating components. We assume the time averaged velocity is 
zero, and there is a very small mean stratification with squared 
buoyancy frequency -w^. Because the convective region is nearly 
adiabatic, aic ^ A'o- We treat the fluctuating components of the ve- 
locity and entropy in the convective region as source terms in the 
wave equation. In practice, we only include source terms due to the 
Reynolds stress in our analysis; source terms due to entropy fluc- 
tuations are of the same size or smaller than the Reynolds stress 
terms (P81, GK90). 

With these assumptions, the equation for the evolution of the 
vertical displacement f . is 

" ...... . (1) 



v'^f. + ^o'vi^^ = o. 



in the radiative region, and 



oz 



(2) 



power spectra. 



in the convective region. We take Vj^ = d^e^, + dyCy to be the hori- 
zontal part of the gradient operator (perpendicular to gravity), and 
S to be the source term due to the Reynolds stress F = V • (u^u^) 
associated with convective motions u^. 

We now discuss the wave solutions to the homogeneous equa- 
tions, i.e., taking 5=0. Because the equations are autonomous in 
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X, y, t, we can Fourier transform in these directions. Thus, we can 
take the solutions to be 

^.(x, y, z, t) = ^-Xz) exp(ifc,-x + - i(jjt), (3) 

and define the horizontal wavenumber kx_ = + k^,, i.e., the 
wavenumber perpendicular to gravity. Throughout this paper we 
will assume A^o » oi. The solutions to eqns. lll2l are 

= Si cos(Nokj^z/co) + Bi sin(A'o^±z/w), (4) 
= Ci exp(-fcxz) + Ci exp(fcj^z), (5) 



respectively, where we have assumed ^o)^ + oil ~ oj. The horizon- 
tal displacement and pressure perturbation Sp are related to 
by 



dp ~ ipo{Now/kj^)^-, 
in the radiative region, and 

6p ~ po(co^/kj^-. 



(6) 
(7) 

(8) 
(9) 



in the convective region. The background density is po, which is 
constant to lowest order in the Boussinesq approximation. 

To solve for the coefficients in eqns.|4]&|5]and the eigenvalues 
o), we must impose four boundary conditions and a normalization 
condition (the latter is discussed in Section [33] >. Two of the bound- 
ary conditions are on the behavior of f . at z = ±L. The physical 
solution requires that f - = at the top and bottom boundaries. The 
other two boundary conditions are set at the radiative-convective 
interface, z = Zi- These depend on the nature of the boundary be- 
tween the radiative and convective regions, and determine which 
0) satisfy the eigenvalue problem. Assume that A'" varies from 
to -a>l in a thin layer with height d, as illustrated in Figure [T] If 
there is a sharp transition between the radiative and convective re- 
gions, i.e., (k±No/cii)d <c 1, we can make the approximation that 
A^^ is discontinuous at z,, which we take to be at z = 0. However, 
if varies slowly, i.e., (k^No/(jL))d » 1, then interesting behavior 
can take place in the transition region. As we discuss in Section[6] 
we expect the most efficiently excited waves in the Sun to fall un- 
der this latter regime. We will discuss both the discontinuous and 
smooth A'- limits below. 



3 WAVE GENERATION BY TURBULENT CONVECTION 

Because the wave generation and wave propagation regions are dis- 
tinct, we use a Green's function (or equivalently, variation of pa- 
rameters), as in P81. Once we have a Green's function G(z, f; (,t), 
we can write the vertical displacement in the radiative region as 



^cad = j drj di;G(z, 



t-(,T)S{x,y,(,T), 



(10) 



where we assume that f,jad is zero at f — » -oo. The Green's func- 
tion depends on whether can be modeled as discontinuous or 
smooth at the radiative-convective boundary. In Section [3T1 we cal- 
culate the Green's function assuming A'^ is discontinuous (as was 
assumed in GLS91 and GK90) and then in Section [3l2l we treat the 
smooth A'^ case. As we shall argue, the latter is more appropriate 
for the low frequency waves which dominate the IGW flux. In Ap- 
pendix |C] we show that the Green's function method is formally 
equivalent to GK90's method of expanding ^- into normal modes 
to solve eqns.[T]|2] 



3.1 Green's Function for Discontinuous A' 

To calculate the Green's function, we need two linearly indepen- 
dent solutions, one which satisfies ^j(-L) = 0, and one that satisfies 
<f;(-l-L) = 0. The boundary conditions at z,, which we take to be at 
z = 0, when is discontinuous, are that <f, and 5p are continuous 
at z = 0. The first solution, which we call rfi, satisfies the boundary 
condition at z = +L: 



''^ " \ Bi exp(-^xz) 



z < 0, 
z>0. 



(11) 



Here we use superscript d to denote the eigenfunction when A'^ 
is discontinuous at the interface. In Section lT2l we will use super- 
script s to denote quantities when A'^ is smooth at z = 0. The sec- 
ond linearly independent solution, which we call ^'!, satisfies the 
boundary condition at z = -L: 



exp(-^xz)) 



z<0, 
z>0. 



(12) 



The eigenvalues a> must satisfy sm{Nokj_L/ a>) = 0. Later we will 
project the total vertical displacement in the radiative zone onto 
the basis {^-]^. The vertical displacement in the radiative zone is 
approximately orthogonal to [r]-]^ in the radiative zone. Thus, it 
is important that our second linearly independent solution is also 
approximately orthogonal to as is the case for egns.ll 1I&I12I 
The general expression for the Green's function, assuming z < 

is 



G(z,i;(,T) 



jda.' 



6(f(co'))^,(z;aj')rj,(^;co') 
w'2 W(0 



exp(-!(jj'(f - t)), 
(13) 

where we label the eigenfunctions with their frequency a>', 5 de- 
notes the Dirac delta function, and W(() denotes the Wronskian 
of ^, and T]-. f(co') is a function which is zero if and only if at' 
is an eigenvalue. For the discontinuous case, we have /(w') = 
sm{Nokj^L/LL)'). We thus can simplify egn.lOlto 



G(z,nf,r) = J] 



1 ^,(z;oj')r,-S(;w') 



Nok^L 



mo 



e\p(-ici>'(t - t)), (14) 



where the sum is over the eigenvalues to'. For the discontinuous A'- 
problem, assuming z < and ^ > 0, the Green's function is 



G\zj;(,r) = j] 



1 



N^klL B2 



^-(z; oj')exp(-k^^-ioj'(t-T)). (15) 



3.2 Green's Function for Smooth A'^ 

If A'^ varies smoothly from A'^ to -lj^, then a WKB type approx- 
imation can be used, provided that NQk±d/a) » 1. Our motivation 
for studying this limit is that the largest scale eddies in stars sat- 
isfy Nokj^d/ci) » 1 (see Section[6}. We would like to develop an 
approximate solution which is valid within the transition region, 
allowing us to connect the solution in the radiative region (eqn.|4ll 
to the solution in the convective region (eqn.|5}. 

The solution in the transition region depends on the form of 
A^^(z) near the radiative-convective interface. In Appendix [B] we 
consider a piecewise linear A'^ profile. This is the most abrupt con- 
tinuous transition possible. We find that the Green's function for a 
piecewise linear A'- is larger by a factor of (A'o/cj_rf/<y)''* than the 
Green's function for discontinuous A'^, implying that wave excita- 
tion is more efficient by a factor of {Nok±d/aj)''^ . 

Here we will consider the case of a tanh profile. We believe 
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a tanh profile is a good approximation to A'- near the radiative- 
convective interface. An eigenmode with frequency oj transitions 
from oscillatory behavior to exponential behavior at a point z, 
(where A'^ = w^), which is lower than the the radiative-convective 
interface, z, (where A'^ = 0). For a tanh profile, z, does not change 
very much as uj changes; although it is smooth, it is not too smooth. 

One might be tempted to appeal to WKB analysis to solve 
for the eigenfunction on either side of the interface, and then match 
across the interface by expanding A'^ to linear order near the wave's 
turning point (as is standard in, e.g., quantum mechanics). Roughly, 
a WKB solution is valid if the local wavelength of the eigenfunc- 
tion is small compared to the scale on which the wavenumber of the 
eigenfunction varies, which for us is d. For smooth A'^ we have as- 
sumed Nok^d/o) » 1, so the WKB solution in the radiative zone is 
valid. However, the WKB solution might break down near the con- 
vection zone if k±d <k 1. Indeed, for the tanh profile, the eigenfunc- 
tion is poorly approximated by the WKB solution when k^d <K 1. 
Moreover, because d/H < 1 and IGWs with k^H ~ 1 dominate 
the wave flux (see Section [331 H here is the pressure scale height 
which we assume is the largest scale of the turbulence), the WKB 
solution fails for the most efficiently excited IGWs. Instead, we 
need to develop a different method to solve for the eigenfunctions. 
The details of this calculation are given in Appendix [A] 

We assume N^(z) is given by 

nHz) = ^° ^ (tanh + l) - oj^. (16) 

In Appendix IaI we derive approximate forms for two independent 
eigenfunctions, and show that there is excellent agreement between 
the numerical solutions to the eigenvalue problem and our asymp- 
totic Bessel function solutions. We are interested in the behavior 
of the eigenfimctions near the radiative-convective interface z, . The 
interface is at 

\ d/ Ni + (ul 



The two independent solutions are 
f Bicos(Nok_i^z/co + n/4) 



z> -d 



B2 s'm(Nok^z/w + 7T/4) z>-d 
5.(^)"^(fr>'.[^exp(-^)] z<d 



(18) 



(19) 



wheretZ)^ = 01^ + a;- and ranges between V2 and 1 for a; > 01^., 
and d = cokj^d/cj. In eqns. [T8]&[T9]we have dropped several fac- 
tors of order unity from the equations derived in Appendix [A] The 
eigenvalues for this problem are the frequencies a) which satisfy 
sini-Nok^L/dj + n/4) = 0. 

Given eqns.[T8]&[T9l the Green's function, for z < and ^ > 

is 



X J f,(z; co') exp(-fcx<r - i<^'U - t)), 
where we introduce the shorthand 

/ aJckj_d\ 

J = Ja 



(20) 



(21) 



Using series expansions from lAbramowitz & StegunI l ll972h . we 
can approximate 



1 



if k^d « 1, 



(k^d)-'" exp(-kj_d) lik^d » 1. 



(22) 



Note that although the smooth and discontinuous Green's functions 
are equal to one another when Nokj^d/a> ~ 1, the Green's func- 
tion in eqn. [20] is no longer valid when Nok^d/o) « 1 (see Ap- 
pendix lAlT l. Instead, eqn.[T5]must be used in this limit. 

3.3 Amplitude Equation 

Now that we have the Green's function, we can calculate mode 
excitation. First, we will expand ^,jad (in eqn. llQt into eigenmodes 
f;,rad(z; w). We usc the subscript rad to denote the z < z, part of the 
eigenfunctions (eans.[T2lll9t. We write 

4'c.iad = — ^ y, A{t; co') ^jjad(z; (^') exp(;l.,x -1- *,,y - ico't), (23) 

where the ^;jad are the z < z, part of the eigermiodes Using this 
representation in eqn.[TO] we take the inner product with f ,_rad(zi '^), 
multiply by exp(-ik^x - ikyj + iait), and integrate over dxdy to find 



A(f ; <jj) ■■ 



X S (x, y, t^, t) exp(-!fcj.x - ikyy + iujT). (24) 

This procedure is discussed more thoroughly in Appendix Icl 

At this point we must pick a normalization condition for our 
eigenfunctions. The energy in the perturbation is 



d^xp^ 



^,ad(z) 



/ 

Y^Yj ^M^*(^') [^^' \^ dzpof,ad(z.'^)-Cdfc^')j(25) 

We want to identify l^('^)P with the energy, so our normaliza- 
tion condition is 



Ww' ^ dz Pof,ad(z; ^) • fr*ad(z; ^') = 



(26) 



where 5 is the Kronecker delta. Using the eigenfunctions (egns.fTTI 
ll2l|18|[T9t and the polarization relation (eqn.|6), the normalization 
condition implies 

B'-~B\~bI~ -i— , (27) 

for both smooth and discontinuous A'^. Using this normalization in 
eqn.|24l the amplitude equations are 



dxdy exp(-!fc,x — ikyy + iwr) 



Nnkl V i J J, 



d( exp(-fcj_f ) 5 (x, V, (, t). 



(28) 



dxdy exp(—ikj,x — iky.y + iair) 



(3)' x 



d^exp(-k^OS(x,y,(,T). (29) 



3.4 Model of Turbulent Convection 

To make further progress, we need to specify the source term 5 . 
We assume that the convective turbulence is comprised of a large 
number of incoherent eddies, estimate the wave generation due to 
a single eddy in isolation, and then find the total wave genera- 
tion by summing over all eddies. We model the statistical prop- 
erties of stellar convection using Kolmogorov turbulence (see, e.g.. 
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iGoldreich &Keelevl[l977h : the convective velocity on the outer- 
scale H is and the associated convective turnover frequency is 
Li)^ ~ u^/H. The convective energy flux is Fconv ~ Po"c- scales h 
sufficiently small compared to H, the turbulent power-spectrum is 
given by the Kolmogorov scaling: 

ui, ^ u, (h/H)"^ ^ u,(coJiD,)-^'^ (30) 

where we have used the fact that smaller eddies have higher fre- 
quencies, i.e., shorter turnover times, with a)^ ^ uiJh oc h^^'^ and 
thus h oc oj/'^'. A given convective eddy characterized by its fre- 
quency 0)^ can excite waves having frequencies o) and horizontal 
wavenumbers that satisfy 



(D<ijj, and fc, < k"' 



(31) 



3.5 Energy Generation Rates and IGW Fluxes 

We can now estimate the energy generation due to a single eddy 
with size h and turnover frequency aie. The source term contains 
three spatial derivatives which we can integrate by parts. The con- 
tribution due to the source term is 



5 ~ t] ui 



(32) 



Assuming the eddy has volume and lasts for a time oj^ ' , we can 
estimate the change in the amplitude due to a single eddy 



AA''(w) 



-k±h ui,. 



3nergy 



(w). 



The total energy generation rate due to all eddies is then 

,.3;,3 



(33) 
(34) 

(35) 
(36) 



The factor of J{k]} jh^ in eqn.[35]counts the number of eddies with 
size h which excite IGWs with frequency cj. We have assumed ex- 
citation happens in a region with thickness dz ~ k]} (because the 
IGW eigenfunction decreases in the convection zone over a char- 
acteristic lengthscale ~ fe^'). Because of the random phases of the 
convective eddies, the excitations due to different eddies are as- 
sumed to be uncorrelated, and the energy increases only linearly 
with the number of eddies. 
The IGW flux is given by 

dF'' E''{oj) / , Afo\ 



d log (jj d log kx_ y\. 



-13/2 



(37) 



-15/2 



2kj_d 



(38) 



d log oj d log k± 

where M = oJcINt, is the convective Mach number. The term in 
parentheses in the first equality of eqn.[37]is the density of states. 
There are J{k\ modes in the horizontal direction, and Lk^^Noloj 
modes in the vertical direction, with wavenumber ~ kj^ and fre- 
quency ~ 0), which each contribute a flux E(a>)/^. Recall that 
eqns. [37] [38] only apply for ai > 01^ and kj^ < ^""(61;) ~ 
H-'((u/oj,f'-. 

Now consider how the flux varies with k^. In the case of 



smooth N^, the flux decreases exponentially for kj^d » 1 because 
of the J- term in eqn. [38](see eqn.l22t. The dominant contribution 
to the flux is from kj^d < 1, so for the rest of this section, we will 
assume k±d < 1, and thus J ~ I. Assuming that d/H < 1, and 
integrating over k^, we find 

\-l/2 



~ poulM 



dp" 
rflogw 

dF' { Po"'(^) 



(39) 

CJ < o)AH/df'\ 
rfl^ " \ POM? {^y'"' (f CO > coAHIdfl\ ^^""^ 
Finally, we find that the total flux is 



F'' ~ pQulM- F,,^,M, 



F ~ Pqu^ 



H 



H 



(41) 
(42) 



This estimate predicts, for a smooth A'- profile, an IGW flux 
only slightly smaller than the convective flux. However, as we now 
show, energy-bearing waves in the smooth A'^ case will undergo 
vigorous wave-breaking as soon as they enter the radiative zone, 
and thus cannot be considered waves. Instead, these large scale 
waves can be thought of a manifestation of penetrative convection 
in this calculation. 

To quantify this argument, we calculate the typical size of the 
perturbations in the radiative zone using 

dF , 

— — — ; po{cj^±)-Ug,z, (43) 

d log cad log k^ 

where Ug- ~ o/k- ~ c? UNok^) is the vertical group velocity, and 
we have assumed s> ^- (eqn.[6ll. From this, we find 

-21/4 



H^ik.Hf" 
No 



and 



kJl ~ M-"'(k^H) 



-21/4 



-23/4 



H 



1/2 



(44) 



(45) 



(46) 



(47) 



Recall that the condition for wave breaking is k-^- ~ 1. For the 
case of discontinuous A'^, the most efficiently excited waves are 
marginally susceptible to wave breaking. However, for the case of 
smooth A'^, the most efficiently excited waves will immediately 
break once they enter the radiative region. Thus, these perturba- 
tions are not actually waves, as they do not stay coherent. 

The only waves that successfully propagate in the radiative 
zone have k-^- < 1. Thus, to find the IGW flux for smooth A'^, 
we must integrate the flux only over the regions of (k±,LL)) space in 
which k,^^ < 1 . This implies 

\-23/4 

(48) 



(HM/d)"'- <{k^Hfi— 
\cJc 



and, as before, a> > o^, kx_ < H^'-io/o^y^-, and k^d < 1. We find 
that the waves that are marginally susceptible to wave breaking, and 
which maximize the flux are at the convective turnover frequency, 
o ~ o^, but have small wave numbers, kj^H ~ (MH/dy^^ . For 
these waves. 



3/8 



(49) 



© ???? RAS, MNRAS OOO.fTlfni 



6 D. Lecoanet & E. Quataert 



This result is only valid if these waves see a smooth A'~ profile, i.e., 



(MH/d)-'"^ » 1 



(50) 



If this condition is satisfied, then the IGW flux is larger than that 
predicted by the discontinuous result by (MH/dy^''^ . Note that if 
d/H ~ M, then the discontinuous and smooth limits give the 
same wave flux. 



4 PRESSURE PERTURBATION BALANCE 

A more heuristic way to derive the IGW flux is to compare the pres- 
sure perturbation on either side of the radiative-convective bound- 
ary. This argument is not sufficiently precise to treat the smooth 
case — hence, we will assume is discontinuous, and thus that 
the pressure perturbation is continuous at the radiative-convective 
interface at z = 0. The pressure perturbation associated with a con- 
vective eddy with a turnover frequency and size h is 

Spcom ~ Povl ~ PouliujJ CDcY\ (51) 

The polarization condition (eqn.|7) relates the pressure perturbation 
in the radiative zone to the vertical displacement. 



Spmd ~ PO- 



(52) 



We assume the convective eddy can only effectively couple to an 
IGW if the frequencies and horizontal wavelengths match, which 
requires 0)^ ~ co and ~ h^' . 

A large number of convective eddies contribute to driving a 
given standing IGW. This is particularly true for k±H » 1 and/or 
ll) » because then small eddies with sizes h <s: H are responsi- 
ble for the driving. The number of eddies contributing to the exci- 
tation of a given standing wave is 



9/1 



(53) 



where we have assumed cj > ojc and that the excitation happens in 
a region with thickness dz ~ k']^ (see also eqn.|35t. Because an in- 
dividual IGW is excited by many uncorrelated eddies, the effective 
pressure fluctuation driving a wave is reduced by a factor of 
relative to that given in eqn. IST] 

When A'" is discontinuous at z = one of the boundary condi- 
tions is that Sp is continuous at z = 0, so that Sp^^^ ~ Sp^om- Using 
eqns. I51I53I we find that the amplitude of a mode with frequency 
~ cl) and wavenumber ~ k, is 



^Pcomk± 



H-k, 



N 



-1/2 



(54) 



However, there are Jlk\ such modes in the domain (we have 
already implicitly summed over the vertical modes in deriving 
eqns. 1511152b . so the typical rms vertical displacement is 



^~ Nodj'- \ N No 



-21/4 



(55) 



the same result as in the inhomogeneous wave equation calculation 
(eqn.l44l(. 



5 COMPARISON WITH PREVIOUS WORK 

In this section, we discuss the relationship between our results and 
previous calculations in the literature. We begin with GK90, who 



only consider the discontinuous A'" case. GK90 solved the fully 
compressible inhomogeneous wave equation by expanding the per- 
turbation in terms of normal modes and then deriving an ampli- 
tude equation. This is equivalent to our Green's function method 
(see Appendix let. Their end result is very similar to our own; for 
k^H \ they find 

\-13/2 

(GK90;eq. 73) (56) 



dF 



--Mpoulik^Hfl — 
d log (J) d log kx_ \ (iJc 

This differs from our result (eqn.l37t by a factor of fcx^^EI We ar- 
rive at a different IGW flux because in the Boussinesq approxi- 
mation d,^, ~ k^^-, whereas for the fully compressible system, 
d-^- ~ ^JH when k^H <c 1. Accounting for both k^H > 1 
and ki_H < 1, the correct scaling of the IGW flux with kj^H is 
F ~ {k^HYii + k±H). Taking this into account in our smooth 
A'^ calculation, the flux of IGWs which do not break changes to 
F'' ~ F^anvM'^'^^(d/Hf'-\ which is very close to the flux in 
eqn.|49] 

Because GK90 solve the fully compressible equations, they 
include multiple scale heights in their convection zone. They find 
that the most efficient excitation of waves with frequency w is at the 
height where the turnover frequency of the energy bearing eddies 
is about equal to o). This effect would be straightforward to include 
in our model — one would need to derive a Green's function based 
on the fully compressible eigenfunctions, and then convolve with a 
vertically varying source term. 

GLS91 use a pressure balance argument to study the discon- 
tinuous A'' case. Their power spectrum agrees with eqn. 1371 when 
Ll) ~ cl)c and k±H ~ 1, but not at higher frequencies or wave 
numbers. They assume that the pressure perturbation in the con- 
vection zone equals the pressure perturbation in the radiative zone. 
They take (5pconv ~ pou\, and (5p,ad ~ Po(wf±)^- This expression 
for the pressure perturbation in the radiative zone does not satisfy 
the polarization condition Spai ~ po{NoLL)/k±^)^- (eqn. |7j, unless 
fx ~ ^l'- Because many eddies contribute to the excitation of a 
single IGW mode, GLS91 also decrease their IGW amplitude by 
a factor of 1 / V^. However, they only account for the incoher- 
ent sum of small eddies at the interface producing perturbations on 
large spatial scales. This gives A/gls91 ~ (k±hy^, where k^h « 1. 
In our analysis, we include eddies which are a distance k^^ above 
the interface, and we take into account that IGWs excited in dif- 
ferent parts of the domain incoherently interfere with each other as 
they propagate in the radiative zone. These additional effects yield 
N ~Jlk-J/h\ 

P81 uses two different techniques to calculate the IGW flux. 
The first uses a pressure balance argument. Press uses that 6p^om ~ 
Poi'l, and that Spnd ~ (PoNoCL)/kj_)^-, and that these pressure per- 
turbations are about equal at the interface. Throughout his analysis. 
Press assumes k^' ~ h. Thus, Press finds 



{kjif- 1 
^'~~ir^~Y' (P81;eq.75) 



(57) 



This is the same result given by GLS91, and is consistent with our 
own assuming kj^h ~ 1. This is because ^k\l N ~ I when kji ~ 1. 



Although our final resufts are similar, there are some ambiguities in 
GK90's derivation. In deriving their eqn. 48 from their eqn. 45, GK90 ap- 
pear to assume that the Sp are orthogonaf under the weighting function c"^ 
and that j dzpoc^^lSpl^ ~ 1. Both of these are true for sound waves, the 
main focus of their paper However, for IGWs, j dzpoc^'^\Sp\^ ~ A1^, and 
the Sp are only orthogonaf under the weighting function 1 (see AppendixICl 
for further discussion on orthogonaiity). 
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Press also derives this result more rigorously using the method 
of variation of parameters, which is equivalent to using a Green's 
function. In addition, Press considers the case in which N- is con- 
tinuous at the interface. He only treats this case in the limit in which 
o) ~ No, and finds 

(P81;eq. 88) (58) 

k 

the same result as eqn. |57] However, note that if tii ~ A'o, then 
Nok±d/a) « 1, and the smooth result cannot be used (these waves 
see the interface as discontinuous). In addition. Press's use of stan- 
dard WKB matching to treat the smooth A'^ profile is generally not 
applicable (see AppendixlAt. 

Finally, we consider the work of B09. In their paper, Belkacem 
et al. numerically calculate the eigenfunctions for a solar structure 
model, use a convection simulation to specify the source term, and 
solve an amplitude equation in the same way as GK90. It is un- 
clear whether the A'' profile in their solar structure model has a 
smooth transition between the radiative and convection zones — if 
their A'^ profile is discontinuous (Section [37Tt or has an abrupt tran- 
sition (Appendix |B](, they will derive different eigenfunctions than 
for a smooth tanh profile (Appendix [Ajl. These eigenfunctions will 
produce a smaller flux (eqns. [41] IB2U than we predict for a more 
realistic radiative-convective transition (eqn.|49t. 

Another key difference is that Belkacem et al. use an eddy- 
time correlation function Xki^^), which in the notation of this pa- 
per can be written as xi'^' '^e)- This function describes how effi- 
ciently an eddy with size l/k and turn-over frequency oj^ = u^k 
excites a wave with frequency oi. Our analysis implicitly assumes 
X(co',ciJe) ~ exp(-(jj^/(jj^). This Gaussian eddy-time correlation 
function implies that eddies with turn-over frequencies oj^ only 
excite waves with frequencies oj. However, the turbulence in the 
convection simulation in B09 is not well described by a Gaussian 
eddy-time correlation function. Instead, Belkacem et al. find that a 
Lorentzian distribution, ;^'(a;; oi^) ~ (1 + 2{ci) I ojef'y^ , is more ac- 
curate. This indicates that waves with frequency u) can be excited 
by a broad range of eddies. In general, this makes wave excitation 
more efficient. It would be straightforward to generalize our results 
to this Lorentzian expression for;^f(w; u),,). 



6 DISCUSSION & CONCLUSIONS 

In this paper we have calculated the excitation of internal grav- 
ity waves (IGW) by turbulent convection, motivated by the appli- 
cation to stellar convection. We assume that the source term ex- 
citing the IGWs can be modeled by Reynolds stresses associated 
with uncorrelated eddies in a Kolmogorov turbulent cascade. Our 
main results are the IGW fluxes, eqns. |4T]&|49l In particular, we 
predict a larger wave flux than previous calculations for low fre- 
quency waves which satisfy Nok±d/<jj » 1, where A'o is the buoy- 
ancy frequency in the radiative zone, kj^ and to are the horizontal 
wavenumber and frequency of the IGW, respectively, and d is the 
thickness of the transition region between the radiative and convec- 
tion zones. We also reconcile somewhat disparate claims in the lit- 
erature by showing that different methods, such as pressure balance 
arguments and solving the inhomogeneous wave equation, predict 
the same IGW power spectrum when using the same assumptions 
(Sectionllll. 

An IGW with frequency ai sees the transition between the ra- 
diative and convection zones as discontinuous if Nok^d/co <c I. 
In this case, the total flux is F'' ~ F^om M, as derived in past 
work, where F^onv is the convective flux and M is the convective 



Mach number. The most efficiently excited waves have frequencies 
0) ~ the eddy turn-over frequency of the largest turbulent ed- 
dies, and fej^ ~ the inverse of the pressure scale height. These 
most efficiently excited waves are marginally susceptible to wave 
breaking when they enter the radiative region. 

If, however, the transition between radiative and convective 
regions is smooth (i.e., Nok±d/a> » 1), the problem becomes more 
complicated. The IGW flux depends on the structure of the buoy- 
ancy frequency A'^(z) near the transition between the radiative and 
convective regions. We parameterize the transition using a tanh pro- 
file, which is a smooth transition that we believe is a physically rea- 
sonable model. The wave excitation is more efficient when is 
smooth because the IGW eigenfunctions change amplitude rapidly 
near the interface (as originally discussed by P81). The total IGW 
flux is F" ~ F^-a^^/id/H) » F''. Again, the most efficiently ex- 
cited waves have frequencies ai ~ a>c and k± ~ H \ However, 
these waves are extremely prone to wave breaking, as k-^- » 1 in 
the radiative region (e.g., P81). These "waves" will thus break im- 
mediately when they enter the radiative region. We interpret these 
perturbations as a form of penetrative convection — they do not stay 
coherent long enough to be truly considered waves. The flux of 
IGWs that are marginally susceptible to wave breaking (i.e., have 
k-^- ~ I) is F ~ F^„,„M^'^{d/Hf'^. This is larger than the discon- 
tinuous A'- flux by (MH/dy^^^ . The flux is dominated by waves 
with frequencies co ~ aic and wave numbers kj_ ~ H-\MH/dy'^. 



In the Sun, air 



IQ-'No, so M 



10-' (e.g., 



iBrown et alj |2012|). and d is estimated to be ~ O.IH 
I Christensen-Dalsgaard et al.l20I ih . IGWs produced by the energy 
bearing eddies have Nokj_d/ci) ~ 10^, and thus the transition region 
must be treated as smooth. This suggests that the IGW flux in the 
Sun is 



F~F^„^M"[j^ 



3/8 



5 X 10"^ F. 



(Sun) (59) 



about five times larger than the flux in the discontinuous A'" case. 
The flux is dominated by waves with frequencies ~ oj^., and k^ ~ 
0.5H-'. 

The increase in wave flux due to a smooth radiative-convective 
interface is only for waves with Nok±d/tL) » 1, i.e., for low fre- 
quency waves. For certain applications (e.g., helioseismology), the 
flux of low frequency waves is unimportant. In particular, low fre- 
quency g-modes in the Sun and massive stars are strongly damped 
by radiative diffusion and are unlikely to be seen at the surface. 
Thus, the increase in wave flux we predict for low frequency waves 
does not change the expected amplitudes of potentially observable 
g-modes in main sequence stars. 

However, low frequency waves are important for the angular 
momentum transport, mixing, and/or mass loss due to IGWs ex- 
cited by stellar convection. For example, a larger IGW flux may 
increase the predicte d mass loss in the final stages of the life of a 
massive star (Ouataert _& ShioddllOl^ and in Type la supernova 
progenitors tPirqi201 ih . This will be studied in detail in future 
work. 

In order to make a more accurate prediction of the wave flux 
and spectrum, one would need to use a stellar structure model with 
a realistic radiative-convective interface and a better representation 
of the convective turbulence, as in B09. Our results highlight the 
importance of adequately resolving the smooth transition between 
the radiative and convective regions in such calculations (which is 
not necessarily the case in many stellar structure models). A dis- 
continuous or abrupt transition will give a qualitatively different 
IGW flux than a smooth transition. 



© ???? RAS, MNRAS 0Q0.[TlfT2l 



8 D. Lecoanet & E. Quataert 



In this paper we only consider IGWs excited by the Reynolds 
stress associated with turbulent convective motions. However, if A'- 
is smooth, we find that the most efficiently excited modes are very 
susceptible to wave breaking. Wave breaking produces turbulence 
and can lead to additional IGW generation (Fritts 2009). Put an- 
other way, penetrative convection can itself generate IGWs. When 
is smooth, the flux in modes which are unstable to breaking 
is a significant fraction of F^onv; thus the breaking process has the 
potential to excite a non-negligible flux of IGWs. 

Perhaps the most promising way to test the results of this paper 
is through comparison with direct numerical simulations of a radia- 
tive zone adjacent to a c onvection zone (e.g.. [Rogers & Glatzmaied 
I2OO5I : Isrun et alj I2OI ih . Although such simulations typically re- 
quire artificially high conduction in the radiative zone, and it is 
unclear how to best identify IGWs dPintrans et alJl200S . this is 
probably the simplest system in which one can quantify the IGW 
flux generated by convection. We hope that analysis of such simu- 
lations can provide a quantitative test of the theory derived in this 
paper in the near future. 
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APPENDIX A: TANH PROFILE EIGENFUNCTIONS 



We will derive the eigenfunctions for the equation 
where 



(Al) 



^(tanh(-i)+ l)-a;2. (A2) 



The transition between oscillatory behavior and exponential behav- 
ior (where N^(z) = ar) is at z, given by 



■ exp 



H)- 



(A3) 



The eigenfunction in the radiative zone is well approximated by the 
WKB solution. 



^, = Bi{NokJa>f^- fc,(z)-"2cosl I dzk,{z) + nlA 



+B2 iNok^lojY'^ k-(zr"^ sin dzLiz) + n/4j , (A4) 



where we define the vertical wavenumber to be 
liiz) = ki [nH,z)Io? - 1) . 



(A5) 



Near Zi, the WKB solution in the radiative region diverges. We 
wish to derive a new set of functions which closely approximate 
the eigenfunctions for z> Zi- 

In many problems, the WKB solutions near a turning point 



can be asymptotically matched onto Airy functions, which provide 
a connection between exponentially decaying and oscillatory WKB 
solutions. However, we cannot use this approach when kx_d < 1; in 
this parameter regime k^iz) cannot be well approximated as linear 
near z,. Instead, we will show that when kx^d < 1 the eigenfunc- 
tions can be well approximated in terms of Bessel functions. Fur- 
thermore, these Bessel function solutions are also a good approxi- 
mation when k^d > 1. 

To show this, first note that if exp(-2z/d) « 1, we can ap- 
proximate 



N\z) ^ {nI + aj^)sxp{-^]- col 



(A6) 



The solutions to the wave equation (eqn. lAlt for this approximate 
N^(z) function are 



k.d— 



+C2Y, 



ci)k± d/d. 



k^d 



4 



^0 + <^c 



■ exp 



■(-i) 
(-^) 



(A7) 



where J and Y are the Bessel functions of the first and second kind, 
respectively. We have also defined or = ar + tol, where aj/a) ranges 
between V2 and 1 . These approximate the solution for large posi- 
tive z- We can asymptotically match the Bessel functions onto the 
WKB solution in the radiative zone (eqn. lA4t . We will make use of 
the following asymptotic forms for J and Y: 



1 (xY 



y„(x) ~ - 



r(a -1- 1) 

T(a) 12 



provided that < x <s: ya+l, and 



Ja(.x) ~ a/ — cos 
nx 



V r ^ ■ I a:7T 7T\ 

Y,Ax)~ -»/ — sm X - — - - , 
y nx \ 2 4/ 



(A8) 
(A9) 

(AlO) 
(All) 



provided that x » lo-^ + 1/4|. 

We must consider two regimes, depending on the size of kj^d. 
First consider k^d <g: 1. We can use the asymptotic formula for 
large arguments provided that 



fc.rf— ^exp(--)»-. 



(A 12) 



This constraint can be satisfied simultaneously with exp(-2z/d) <c 
1, implying that the asymptotic form of the Bessel functions are 
good approximations to the eigenfunctions. If we approximate A;^(z) 
as 



J.2 + 



exp - 



2z 



(A13) 



we can approximate eqn. IA7l bv 

~ 1.1. 




d(^(z)) 
-rf(^(z))" 



1/2 nGik^d n\ 



n<j)k,d n\ 

+ + - 

2w A] 



(A14) 
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This matches onto the WKB solution in the radiative region since 
kj^d is small. The amplitudes are 



C? — —Bi 



nNokj_d 

nNok^d 
2cj 



1/2 



(A15) 



(A16) 



Now assume k±d » 1. In this case, the asymptotic form 
of the Bessel functions for small argument is only valid when 
exp(-z/rf) » 1, i.e., for positions where the Bessel functions them- 
selves are not a good approximation to the eigenfunctions (N^(z) 
cannot be simplified as in eqn. IA6l if exp(-z/rf) » 1). However, 
in this limit we can use the WKB approximation in the convective 
region, and connect the two WKB solutions with Airy functions. 
Thus, in the convective region, we have 



~ (BJ2)(NokJcjy'\(zr"'~expl^- J~ dz'\k,(z' 

+B2(Nak^/ajy'%{zr"^ exp |+ dz'\k,(z') 
For z much larger than z„ this becomes 



(A17) 



(z - z,)k^^aj/a)) 



+S2(-^) Fj exp(+(z-z,)k^<Ij/oj), 



(A 18) 



For z much larger than z,, the Bessel functions are a good ap- 
proximation to the eigenfunction. In the limit of large z, the Bessel 
functions become 

-1 '^'^Pl 



-C2 



' 1 

^enkx^ 

4 



(z - Zt)k^(jJ 



12 I 2^^ibk^dluj+\l2 

eu) j 



exp 1+ 



(z - z,)k^(l) 



(A19) 



enki^dj 

Thus, the Bessel function solution matches onto the WKB solution 
in the convective region when 



Ci =Bi 
C2 — —Bi 



2aj 



(A20) 



(A21) 



Using eans. lA15llA16l & IA20IIA2n we can approximate ^- by 



/ nNok^d 
2oj 



lnNok±d 
2aj 



1/2 



V^o + ^? / z\ 

k^d exp - — 

o) \ d! 



k±d exp -- 

(1) \ dl 



(A22) 



where we have defined d = ajk^d/o). This will be a good approxi- 
mation for as long as exp(-z/rf) <c 1. 

For the purposes of determining the convective excitation of 
IGWs, we are interested in evaluating ^. between z, and Zj + 1/fcx, 
where zi is the location of the interface between the radiative and 
convective regions. Since 

-^~exp(-2|), (A23) 

the argument of the Bessel functions varies from <jjck±d/a> to 



exp(-(k±dy^)a)ck±d/a>. Within this range, the Bessel functions 
change by about a factor of e. It is thus within the accuracy of 
our calculation to take ^- to be about constant within this range: 
at z = Zi, we have that 

?z~ Bi\ — — — I I — 1 J^k^d/oj (aJck±d/a)) 



2a> 

+B2\ [-] Y^k,d,^(co,k^d/co). 



(A24) 



In evaluating eqn. IA24I we need to calculate Jj,{xa), where x = 
ibkj^d/cL), and a = oj^-Ioj < I/V2. A good set of approxima- 
tions for the Bessel functions for x « 1 and x » 1 is given 
in eqn. l22l of the main text ( based on expansions of J^(xa) from 
lAbramowitz & StegunI (Il972h '). 



Al Numerical Verification 

Here we will present numerical verification of our approximate 
solutions in the above subsection. We numerically integrated the 
homogeneous differential equation (eqn. lAlt with given by 
ean. lA2l in Mathematica using the "ImplicitRungeKutta" method, 
and solved for a physical solution, satisfying ^, — > as z — * 00. We 
pick the right boundary to be a point b deep within the convective 
region, where k^(b) = -k\, specify ^^(b) = I, ^'.(b) = -k^, and 
integrate ^, leftwards into the radiative region. This ensures that 
satisfies the boundary condition z — » +00. We find that our calcula- 
tions are insensitive to the value of b, provided that it is sufficiently 
larger than z,- 

To test the approximations described in the above subsection, 
we calculate the value of the physical eigenfunction at the interface 
between the radiative and convective regions ^-(zd- Because any 
multiple of the eigenfunction is also an eigenfunction, we normal- 
ize by B[ (see ean. lA4l >. which is the amplitude of the oscillations 
deep in the radiative zone. Equation lA24l Dredicts 



2aj 



[■t) J^k^di,.[k^d-^Y (A25) 



Our analysis is only valid if we are in the smooth A'^ limit, i.e., if 

Nok^d/a> » 1. 

In Figure lAll we compare our numerical results to the an- 
alytic predictions. In Figure lAll (top panel) we vary w/A'o for 
two different values of k_id. The numerical solutions agree with 
our prediction when Nok±d/cL) » 1. In the opposite limit, when 
Nok^d/o) « 1, we can treat A'^ as discontinuous, so ^- is contin- 
uous across the interface, and ^-(zi)/B[ = 1, as is the case for the 
lower curve in Figure |AT] (top panel). In Figure |AT] (bottom panel) 
we vary kj^d, fixing lj/Nq = 0.01, for two values of oj^INq. In this 
case, we have Nok±d/a) = 1 when k±d = 0.01. The normalized 
eigenfunctions approach one as k^^d decreases, and the numerical 
solutions begin to deviate slightly from the analytic prediction near 
k±d = 0.01. These results indicate that our analytic solution for 
if- near z, is accurate provided we are in the smooth A'' limit. The 
numerical solutions also show how the eigenfunctions transition 
between the smooth and discontinuous A'^ limits. 



APPENDIX B: PIECEWISE LINEAR A'^ 

In the limit of smooth A'^, the eigenfunctions. Green's function, and 
IGW flux all depend on the nature of the transition between radia- 
tive and convective regions. In this paper, we focus on the case of a 
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Figure Al. The normalized eigenfunction at the radiative-convective in- 
terface Zj. The symbols denote the numerical solution, and the lines denote 
the analytic prediction, egn. IA241 In the top panel, we vary oj/Nq, fixing 
Wc = <jj. The blue Une and crosses have kj^d = 0.1, and the red line and 
asterisks have kj^d = 0.01. The numerical solution matches the analytic 
prediction for smooth when Nok±d/oj 2> 1, and approaches one (the dis- 
continuous solution) when Nok±d/oj « 1. In the bottom panel, we vary 
k±d, fixing w/No = 0.01 and setting oJc/Nq = 0.01 (blue curve, crosses) 
or Wc/No = 0.002 (red curve, asterisks). Again, there is good agreement 
between the numerical solution and the analytic prediction. 

tanh profile (Appendix|Al(, as we think it is the best simple model of 
this transition region. However, in this appendix, we consider an- 
other analytically tractable transition — a piecewise linear pro- 
file. This is the most abrupt transition possible, and thus not very 
physical. However, we find it instructive to consider a second exam- 
ple of a transition to explicitly demonstrate that the results depend 
on the details of the transition region. In particular, this highlights 
that care must be taken to use a smooth A'- profile when numeri- 
cally solving for wave excitation using a stellar structure model (as 
in, e.g., B09). 

We assume A'^ is given by 



^5 



(A'^ - a,l)/2 - (Nl + 0,1) {zld) 



if z ^ 

if - d/2 < z < 
ifz>dl2. 

(Bl) 



We have that A^(z) = or at the point 



A'^ - lay" - col 



and that N-{z) = at 



(B2) 



(B3) 



The solutions in each region are 

= B,cos(Aoytx(z-^d/2)/aJ)-^B2sin(A'ofci(^-^^//2)/al),(B4) 
for z < -d/2, 

= Ctexp(-kAz-d/2)) + C2exp(kAz-d/2)), (B5) 
for z > d/2, 

= D,Ai{Kl'\z-z,)) + D2Bi{K['\z-z,)), (B6) 
for -d/2<z< d/2, 

where Ai, Bi are the Airy functions of the first and second kind, 
and 



dkliz) 



dz 



klNl + a,l 



(B7) 



We can relate the six coefficients in eqns. IB4IB6I to one another 
using four boundary conditions: ^, and must be continuous at 
z = ±d/2. 

First consider the boundary at z = +d/2. The argument of the 
Airy functions at this boundary is 



(k^d) 



2/3 



Ci/kj^d 



2/3 



(B8) 



This is much smaller than one unless kx_d is extremely large. One 
can check that IGW excitation is exponentially suppressed when 
(x)^k±d/Ng » 1. Thus, we will assume that aPhj^d/Ng <c 1. This 
implies that Ai|rf/2> Bi|rf/2> Ai'|f;/2, Bi'|f;/2 are all of order one, where 
we have introduced the shorthand Ai|; = Ai(A'j''^^(z - Zr)), and sim- 
ilarly for the other functions. To order of magnitude, we have that 



Ci+C2~ £>lAi|rf/2 + 02Bi|rf/2, 



and 



Notice that 



C1-C2 



K 



1/3 



7-(£>iAi'|j/2+£»2Bi'|d/2). 



1 NI + ojI 

kid 



2\l/3 



» 1. 



(B9) 



(BIO) 



(BID 



Now consider the boundary at z = -d/2. The argument of the 
Airy functions at this boundary is 



(k^d)- 



13 "0 



2x1/3 



,.2 



Nok^d 



2/3 



» 1, 



(B12) 



where the last inequality follows from assuming that we are in the 
smooth limit. We thus have 



Nok^d 



.2Nok^d n\ . (2Nok^d n 

£»! cos + - + £)t sin h - 

'3 a; 4/ " \3 a. 4 



implying 



Nok^d 



4/6 



{Di cos{(p) + D2 sin{(p)) , 



(B13) 



(B14) 



where (p = (2/3)(Nok^d/cji) + n/4. Similarly, by comparing on 
either side of z = -d/2 we find 



B2~ 



Nok^d 



i-Di sin(0) + D2 cos(0)) . (B 15) 
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Using these boundary conditions, we find that the physical 
eigenfunction is 

^ r B,cos(^^^^i^) + fl2sin(^^^^) z<-d/2, 
~ I (Mt^y' sxp(-kAz - d/2)) z > d/2, 

(B16) 

where we use superscript I to denote the eigenfunction for the 
piecewise linear profile, and 61 ~ B2 ~ B^. An unphysical 
eigenfunction is 



j z < -d/2, 

Bi sxpi-k^z - d/2)) + §2 expikAz - d/2))) z > d/2, 

where B2 ~ B[ ~ Bi- Note that the constants -61,-62 in rjl and 
-61 , -62 in vary sinusoidally with d (as well as the other parameters 
of the problem). Thus, although for most values of d they are the 
same size, there are specific values of d for which one term is much 
larger than the other. 

The Green's function for z < and <f > is then 

ifNokiy/L\ ^' 
X f'(z; uj') exp(-fcj_^ - icL>'(t - t)). 
This implies that the IGW flux is 



(B18) 



d log o) d log kj 



^p,ulM"Hk.H)'i- 



(k^d) 



1/3 



(B19) 



assuming co > ai,., kj^ < ^"'"(a;) ~ H '(01/01^)^'-, and kj^ 
d^HNo/(jS)^ . The typical rms vertical displacement satisfies 



kS. ~ M-"''(k^HY"'\ — 



-65/12 



(B20) 



We find that the most efficiently excited waves are susceptible to 
wave breaking. The flux of waves marginally susceptible to wave 
breaking is 



F' 



3/22 



(B21) 



Using M ~ 10-3 ^j^^ ^ ^ Q (Section [6l(, we find F' ~ 
1.87FconvA4. This is a very modest increase over the discontinuous 
IGW flux (see Section[33]l. 



APPENDIX C: MODE PROJECTION FORMALISM 
(GK90) 

In GK90, an amplitude equation is derived by projecting the in- 
homogeneous wave equation onto specific modes. We will show 
that their approach gives the same result as our Green's function 
approach, provided that the correct inner product is used. 

First start with the inhomogeneous equation for ^- in the 
Boussinesq approximation 



(CI) 



In the mode projection formalism, we decompose ^, as 

4 = — ^ Mt', (J->')nziz; w') exp(ik.tX + ikyV - icj't), (C2) 



where ri^(z;u)') are the physical solutions satisfying the homoge- 
neous wave equation. Substituting this into the inhomogeneous 
wave equation, multiplying by p(,riZ{z;cj)e,xp{—ik^x — ikyy + icot) 
and integrating over d^xdt, we find 



m-oj)\ 



2k] VI 



/:*/ 



dxdy exp{—ikxX — ikyy + iwr) 



X J d(poS(x,y,(,T)Ti:((;cj). (C3) 
A crucial step in deriving this is using 



dzpod-ri-(z;cij')d-ril(z;cj) = 6^ 



k] 



(C4) 



That is, the j/.(z; oj) are orthogonal with respect to the inner prod- 
uct {a,b) = j dzpod-ad,b' . This follows from our normalization 
equation (eqn.l26b and the polarization conditions (eqn.|6j. 

Although we use ^- as our perturbation variable in this paper, 
GK90 uses 6p. The inhomogeneous wave equation for Sp in the 
Boussinesq approximation is 

V^^_Sp + N'vl6p = S, (C5) 

where S ~ (poor/k^) 5. As above, we can decompose 6p into 
eigenmodes 

6p = — — V' A(f; CL)')Sp(z', co') exp(ikyX + ik^y - iuj't), (C6) 

where 5piz', (jj') are the physical solutions satisfying the homoge- 
neous wave equation. When we put this into the inhomogeneous 
wave equation, multiply by p^dp'iz; oj) expi-ikj^x- ik^y + ia)t), and 
integrate over d^xdt, one might think that 



\A{t;w)\ 



— f 

Ipl yfM J z, 



dxdy exp(—ikxX — ikyy + iwt) 



d(poS (x,y,(,T)6p'{i;;oj). 



(C7) 



2oJNyo^^M. 

Using Sp{(;(jS) ~ {poCiJ^/k^)T]-((;(u) (eqn. |9]l, we see that this es- 
timate of \A(t; ll))\ differs from our estimate using ^- (eqn. IC3t by 
cj^/Nq. This leads to an underestimation of the flux in IGWs by 

The discrepancy is due to using the incorrect inner product. 
Implicit in the derivation of eqn. IC7l is the assumption that the 6p 
are orthogonal under the same inner product as the i.e.. 



dzpod-dpiz; oj')d-6p'{z; cn) = 5^j^> pINq. 



(C8) 



However, one can check that the 6p are not orthogonal with respect 
to this inner product^ Rather, they are orthogonal with respect to 
{a,b} = j dzpQ^ab', i.e.. 



/ 



dzpQ'6p(z;co')Sp'(z;ci>) = 



kY 



(C9) 



3 Using the properties of Hermitian operators, one can show that the Sp 
IGW eigenfunctions of eqn. lCSI are orthogonal under the inner product de- 
fined in eqn. |C8| However, for the mode projection to be well defined, we 
must work in a complete basis, and the IGWs alone do not form a complete 
basis (in the convection zone). Our resolution of this apparent inconsistency 
is to note that the eigenfunctions of the full non-Boussinesq wave equation 
do form a complete basis (this includes sound waves in addition to IGWs). 
Moreover, one can show that the Sp eigenfunctions for the non-Boussinesq 
equations are only orthogonal under the inner product defined in eqn. lC9l 
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Thus, if we integrate the inhomogeneous wave equation twice with 
respect to z, multiply by p^^6p'(z',(jS)e.xp{-ik^x — ikyj + ia>t), and 
integrate over cPxdt, we get 



\A{t-aj)\ 



1 



2aj3 



^ rfT ^ dxdy exp(-ikjX - iky.y + icot) 

X ^ rf^Po'5(x,j,f,T)<5p*(<r;^^). (CIO) 



One can check that this is consistent with the calculation using 

If one uses a Green's function this issue of orthogonality under 
different inner products becomes trivial. Using the expansions in 
Sec. 13. 31 we have 



— = V A{t\ w')f-iad(z; w') exp(iliX + ikyj - iai't) = 



Nok^LWiO 



5 exp(- W(/ - t)), (Cll) 



where z < Zi- Since both the left and right hand sides are in the span 
of (^;,iadLi we can simply use the inner product defined by 



(CI 2) 



Taking (^;,rad(z; <^), ■) of eqn. lCl II multiplying by exp(-ik^x—ikyy+ 
ia>t), and integrating in the horizontal directions, we get 



A(t; oj) ■■ 



X S{x,y,^,T)sxp{-ik^x— iky.y + iair). (C13) 



This is eqn. [24] which can easily be manipulated into eQns. l28l|29l 
using the eigenfunctions. Note that we cannot use such an inner 
product in the mode decomposition formalism because we need 
to calculate terms like ((5/?*(f;aj),5(x, y, <f, r)>, and thus need an 
explicit formula for the inner product in terms of integrals over 

Finally, we will demonstrate that the mode projection 
formalism — when done correctly — and the Green's function for- 
malism give the same result. Specifically, we will show that 
eqns. |C3| and |C13| are equivalent. First note that W{^) is a constant 
for our wave equation. We want to show that 



1 



2kl' 



Nok^LW 

We can evaluate W in the radiative zone, and find 

INok^B^Bi 2ki_ 



W - 



NoCLiLpQ ' 



(C14) 



(C15) 



where we have used eqn.|27] This proves that the two formulations 
are equivalent. 
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